Preprint typeset in JHEP style - HYPER VERSION 



hep-ph /06xxxx 



A Consistent Prescription for the Production Involving 
Massive Quarks in Hadron Collisions 



Borut Paul Kersevan"' 

a Jozef Stefan Institute, Jamova 39, SI- 1000 Ljubljana, Slovenia 
b Faculty of Mathematics and Physics, University of Ljubljana, 
Jadranska 19a, SI-1000 Ljubljana, Slovenia 
E-mail: 



borut .kersevanOijs . si 



Ian Hinchliffe* 

Lawrence Berkeley National Laboratory, Berkeley, CA , 94-720 USA. 



E-mail: I_Hinchliffe@lbl.gov 



Abstract: This paper addresses the issue of production of charm or bottom quarks in as- 
sociation with a high process in hadron hadron collision. These quarks can be produced 
either as part of the hard scattering process or as a remnant from the structure functions. 
The latter sums terms of the type (a s log(pT/m q )) n . If structure functions of charm or 
bottom quarks are used together with a hard process which also allows production of these 
quarks double counting occurs. This paper describes the correct procedure and provides 
two examples of its implimentation in single top and Drell-Yan at the LHC. 



Keywords: |QCD, NLO Computations, Parton Model, Hadronic Colliders 



'Work supported by the Director, Office of Science, Office of High Energy Physics, of the U.S. Depart- 
ment of Energy under Contract DE-AC02-05CH11231. 



Contents 



[l]. Introduction |l| 

[2]. Kinematic Issues ^ 
|2.1| Phase-Space Transformation H 



2.2 Monte-Carlo Generation Steps 



[3L Examples of the Algorithm Implementation 12 
|3.1| Associated Drell-Yan and b-quark Production 12 



3.2 The 't-channel' Single Top Production Process 14 



3.3 The 'tW channel' Single Top Production Process 17 
Conclusion 19 
|A]. Kinematic Relations 19 



1. Introduction 

The production of states by a hard scattering QCD process in hadron collision is described 
by the parton model which separates it into perturbative (calculable) process and a soft 
non-perturbative piece. Consider the production of a top quark pair via the partonic 
process gg — > it. There are many gluons present in the process but the only top quarks 
arise from the hard scatter itself. Top quarks can also be produced singly via the process 
gb — * W~t. In order for this to occur the incoming bottom quark is viewed as a constituent 
(parton) of the incoming hadron. Alternatively one could begin with a two gluon initial 
state and consider the hard process as gg — > tbW~ . These two processes cannot be added as 
the QCD approximation that produces the b parton in the first case is partially accounted 
for by the latter process. A careful separation of "hard" and "soft" components is needed 
so that a consistent result can be obtained. The rest of this paper demonstrates such a 
separation. The rest of this section is concerned with introducing the formalism. Section 
2 shows explicitly how to relate processes with one more (or less) parton in the hard 
scattering. Section 3 presents some explicit examples of the Monte-Carlo implementation 
of this formalism. Finally some conclusions are drawn. 

The Q CD-based parton model is based on the factorization theorems @, H|, [|, ||| 
according to which the squared amplitude for a process A, B — > X can be decomposed into 
a "hard" and "soft" (or alternatively denoted "short" and "long distance") parts: 

\M A B^x\ 2 = fa/A ® H ab^X ® fb/B 
a.b 

= ^f'Y L j -^-fa/A(£,a,fJ-F) fb/B^b, I^f) H ab ^ X ((,a, £b, ■ ■ •), P- 1 ) 
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with a, b labeling the incoming partons which have to be summed over and H(ab — * X) 
denoting the hard ('short time') part of the squared amplitude. The soft contributions 
are absorbed into the parton distribution functions fi/i(£i>MF) with fj,p being the (factor- 
ization) scale at which the two parts were separated 1 . More explicitly, the above theorem 
states that the collinear (mass) singularities have to be isolated/subtracted from the hard 
process amplitudes and reabsorbed into the parton distribution functions 11, BL & [6|, M; 
all other singularities (UV, soft IR) appearing in the perturbative calculation of the hard 
process either cancel or are handled by renormalization procedures. It has to be stressed 
at this point that the renormalization/regularization scheme used in subtracting the UV 
singularities in turn dictates the precise form of the evolution (DGLAP) equations of the 
parton distribution functions Jl], ||| 

1 

J^fi/l(^) = J jP^^F))f 3 ,M^F\ (1.2) 

3 z 

where Pj^i denote the usual DGLAP evolution kernels and I describes either a parton or 
a hadron. 

An elegant way of isolating the mass singularities in perturbative calculations is found 
by observing that the pQCD squared amplitude |A4 a b^x| 2 involving initial state partons 
a, b is subject to the same factorization theorem: 

\M ab ^ x \ 2 = fc/a ® H cd ^ x ® fd/b, (1-3) 

c,d 

with the fj/j representing the distribution function of the parton i inside the parton j. 
The above equation holds to any order in perturbation theory. Consequently, since the 
|-M a 6^x| 2 can be calculated to any order by using the Feynman rules and the prescriptions 
for calculating the fj/j to high orders in a s are also well established [|], ^l], 22] one can use 
the procedure of ]5|, |6|, [7j to extract the H c ^x to the chosen order. 
At zero-th order in a s : 

/$(fl = 3*(e-i) (i.4) 



and hence: 



Subsequently, at first order in a s : 

e (0)m , Al) 



fi/j(o = f^(o+n};io, (i.6) 



and thus at this order: 



In case all partons are considered massless the flux factor in the partonic cross-section expression 



is s = £ a £b(2s) with (2s) being the hadronic flux and the Eq. 1.1 results in the common expression 

CTAB^X = JdCa/dCbfa/A(Ca,MF)fb/B(Cb,MF)cr^ x (s,AlF) 
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where formally the virtuality /x 2 of the particles cord (which is in turn proportional to the 
(pt) 2 of the particles c or d, see Fig. p]) is used as the factorization measure with the limit 
fjLp'. The hard part H^^L^ includes the cases fj? > ^ and the soft one the cases fi 2 < //|. 
Note that in the above equation, the phase space integration over c/ d particles has already 
been performed (resulting in the convolution integral) and thus the phase space for the 
final state particles for the soft part formally involves one less particle. 




Figure 1: The diagrammatic representation of the method applied in isolating the soft (collinear) 
terms in propagators corresponding to virtual particles p c and/or pd when calculating the 
|A^ahLxl 2 - Note that in factorization theorem the particles incoming into the hard part of the 
probability amplitude are considered to be on-shell. 

The above equation can be inverted to give: 

c d 

This can be extended to higher orders in perturbation theory. It should be emphasized that 



the presence of the subtraction terms in the above Eq. 1.8 prevents double counting when 



performing the perturbative cross-section calculation since the collinear effects present in 

the I-MILxI 2 a 
initial hadrons. 



the |-^ib-»xl 2 are removed and re-summed in the parton distribution functions f a M of the 



After the perturbative expansion of H a ^x = Hab-*X + ^a6->X' given by Eq. |L5| and 



Eq. 1.8, is inserted into the cross-section expression (Eq. [O]) one thus obtains the formula: 
\M AB -,x\ 2 = + I-^aLxI 2 - \Mab^x\1 (1-9) 



with the subtraction terms given by: 

\M AB ^ X \l =E/aM®E/i/a®^X®/6/B+E/«M®E^l^®/i>/VB- (1-10) 
a, b c a,b 



For further discussion on the 'double counting' issues it is illuminating to calculate 
the fi/j at scale /if up to the order of a s . One starts by writing down the perturbative 
expansion of the evolution kernels in the DGLAP equations: 
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and observing that the evolution increases the order of f$ by one. This can explicitly 

be seen by inserting the zero-th order fy^ of Eq. |1.4| into the Eq. |1.2| , and increasing the 
factorization scale from the lowest kinematic limit (mass m of the particle i) to the scale 
lip, i.e. integrating over the range [m 2 ,^] and keeping only the terms up to the order of 

/.M5, W )=/iJ?«) + ^^ i (f)m(§). (i.i2) 



The above expression matches the perturbative expansion given by Eq. |l.q with the second 
term identified as the ijj) parton distribution function. 



While the expression of Eq. 1.9 can subsequently be used for estimating the total cross- 
section of a given process one should take further steps when dealing with the estimation of 
the differential cross-sections or (equivalently) Monte-Carlo simulation. In a Monte-Carlo 
simulation the DGLAP parton evolution, re-summed in the parton density functions fo/j 
(c.f. Eq. [E^), is made explicit by evolving the factorization scale from its initial value 
to the lowest kinematic limit (or an imposed cutoff). Identifying the evolving factorization 
scale with the virtuality fx 2 of the incoming particle c the probability that the particle c will 
be un-resolved into a particle a (or equivalently, that a particle a will branch and produce 
the particle c and an additional spectator particle) is given by the Sudakov term (see e.g. 



' ! ' ; ' ^~^^J - P ^ Z) 1^)>- (L13) 

At each evolution step (branching) the number of particles is increased by one and its 
contribution to the differential cross-section in terms of a s is also increased by one. The 
described procedure is commonly known as (initial state) parton showering. 

The (next-to-leading order) subtraction terms of Equation 1.10| thus compensate for 



the first branching in the backward evolution of the incoming partons participating in the 
(leading order) term H^_^ x . Note that in order to match the subtraction terms in Eq. 1.8 



with the first-order matrix element, the fraction £ c (or equivalently £d) of the evolved parton 
is kept constant and the virtuality is decreased corresponding to a 'backward' evolution in 
time from the starting point of virtuality /Xq of H^_^ x . to the virtuality li 2 of H^'_ >x (c.f. 
Figure |l|) . Writing the expressions of Equation l.lCj in differential form in [i 2 one thus gets 
for the first term: 

I— ^fE/f/f/f {w^>^(f) 

a,b,c ' 



H%l x (£ c ,t; b )f b/B (s b , $) + ...} (i.i4) 



and an equivalent expression can be obtained for the second term of Eq. 1.10. Using 



again the Eq. L5, multiplying by the flux factor l/2y A(s, m? A , m 2 B ), given by the Lorentz 
invariant function: 

A(s, m\ , m 2 .) = (s — (mi + m2) 2 )(s — (mi — m2) 2 ) (1.15) 
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and integrating over the final state n-particle phase space denoted by J d&x, one obtains 
the first subtraction term: 

da^(AB^X) _ 
d£ a d(, b dn 2 d£ c d(f) 

|^2lxl 2 (£ c ,6)*&x-c (1.16) 

and equivalently also the second subtraction term by appropriate replacements a — > b and 
c — > d. The two derived equations correspond to the expressions obtained by Chen, Collins 



et al. |9|, |lCj, [11, 12|, derived by the Sudakov exponent expansion. 

In addition to writing the cross-sections in differential form, the integration over an 
angle <j) was introduced, where the angle <p represents the azimuthal angle of the spectator 
particle c or d and is in effect a dummy quantity which nevertheless has to be sampled 
in a Monte-Carlo simulation procedure. The notation d^x-e denotes that the final space 
integral does not contain the spectator particles c or d since they are already accounted 
for in the (d£/£)d^ 2 d</> differential (see e.g. pi]). In contrast the first order matrix element 



of Equation [L8] is integrated over the full phase space X: 
da^{AB -► X) 



dx a dx b 

E / 1 = fa/A( x a, l4) fb/B{?b,$) I \M^ x \ 2 (x a ,X b )d^ X A 1 - 17 ) 

a ,b 2J X(s,m 2 A ,m B )x a x b J 



whereby the outstanding issue is the kinematic translation between n — 1 (X — c) and n (X) 
particle kinematics. Furthermore, one cannot simply equate the variables £ a b and x a b, 



since e.g. in Eq. 1.16 the £ Ci b terms imply that the incoming particles c and b are on 
shell in the matrix element calculation while in Eq. |1.17| in contrast the particle a (as the 
'parent' of particle c) is the on-shell one. 

The mapping of the kinematic quantities between between the expressions of order 



as {i.e. the pQCD derived expression of Eq. 1.17 and the (showering) subtraction terms 



of Eq. 1.16| ) needs a consistent and possibly a formally correct definition. In order to 



achieve this the prescription developed by Collins et al. of how to merge non-perturbative 
(parton-shower) calculations with the leading order perturbative pQCD calculations on the 
level of Monte-Carlo simulations 1C, 11, ^] was implemented, which has been explicitly 



shown to reproduce the NLO MS result for a set of processes. 

Another outstanding issue is that in case of heavy quarks participating as the initial 
state partons the (commonly used) approximation of treating the incoming particles as 
massless can lead to a significant error. This fact, as well as the the solution in terms of 
consistent treatment of the kinematics in terms of light-cone variables, has been demon- 
strated in the ACOT prescriptons of how to consistently introduce the factorization in 
case of non-negligible masses of the colliding partons (e.g. heavy quarks) |], f?J. The ACOT 
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prescription has in this paper been introduced in the formalism of Monte-Carlo simulation 
by modifying the prescription of Collins et al. accordingly. 

In the Monte-Carlo generation steps one thus has to produce two classes of events; one 
class is derived from the leading order process with one branching produced by Sudakov 
parton showering and the second class are events produced from the next-to-leading or- 
der hard process calculation (i.e. the pQCD calculation with the appropriate subtraction 
terms) . 



2. Kinematic Issues 



2.1 Phase-Space Transformation 

The aim of this section is to derive generic ex- 
pressions that transform the kinematics from the 
'hard' to the 'soft' (or showering) n-particle sys- 
tem, i.e. split the 'hard' n-particle phase space 
involving heavy quarks into n — 1 'hard' ©1 'soft' 
particle phase space in order to perform appro- 
priate MC simulation (c.f. Figure |2[). As already 
stated the existing prescriptions deal either with 
massive particles || |7[ on the level of integrated 
cross-sections or with explicit Monte-Carlo algo- 
rithms involving light (~ massless) particles (e.g. 
1 25, |TT| , P~3| , ) while no generic combination 
of the two algorithms is available. 

In order to accommodate the particle masses 
it is convenient to work in light-cone coordinates p^ 



*x(Sn = (k n _ 1 + Pn ) 2 ) 



fcn-i; *£_! = »„-! 




Figure 2: The diagrammatic representa- 
tion of the method applied in translating 
a t-channel (space-like) split going from n 
to n — 1 particle phase space. 

,p r ,p~) where p 1 * 1 = -^(p ± 



p ) and the remaining two coordinates are considered 'transverse' p 



The kinematic 

prescription of relating the hard n-particle kinematics to the soft n — 1 case is as follows: 

1. Incoming hadron A is moving in the z direction and hadron B in the — z direction, 
carrying momenta Pa and Pb with the center-of-mass energy i/s, whereby one can 
neglect the hadron masses at LHC energies. 

2. The incoming partons with momenta p a and pb have the momentum fractions p+ = 
x a P^ and p^ = XbPj^ relative to the parent hadrons and the center-of-mass energy 

(c.f. Fig. D 



3. The split propagator (i.e. particle c) virtuality p\ = (p a — Pc 



gives the fi 2 value from t r 



_ x - m \ = -/i 2 . (c.f. Eq. |Llq) . 



4. The scale correspondence is given by s n _i = /Uq (c.f. Eq. 1.16| ). 

5. All incoming and outgoing particles (partons) are on mass shell. 

6. The splitting parameter of the evolution kernel is z = I s - (already used in Eq. 1.16) 
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7. The rapidity y = ^ In ( J of the subsystem (c.f. Figure ||) is preserved in the 
translation. 



In order to further relate the n and n — 1 particle 

phase space for hard and soft interpretation one can use the recursive t-channel splitting 
relation pi |i~8|]: 



Sn j m\ , TU2 , • • • , Tl n J 

(V^T-m-n) 2 



2~ 



n-1 



4JX(s n ,ml,ml) J Q 



dip* n / dtn_i $ n _i(s„_i,mi,m 2 , . . . ,m n _i), (2.1) 



n-1 
± 



where 99* = (f> of Eq. 1.16 and the limits t n _ x are given by analytic, albeit complex 
expressions [|l(| Using the above relation one can introduce the X — c particle phase 
space split into Equation |1.17| by identifying p E = p n : 

daW(AB -» X) 



dx a dxi,dt n -ids n -idip* 

^ / = fa/A&a, S n -l) fb/B&b, «n-l) X 

a ,b 2J \{s,m 2 A ,m 2 B )x a Xb 



7 4y / A(s n ,m2, 



X-c 



(2.2) 



Using the translation prescriptions introduced above, the remaining issue is the relation 
of the variables £ a ,£o£b with the variables x a , Xb, s n -i- The requirement (2) in the list 
above ensures £ a = x a since both particles in question are incoming partons originating in 
hadron A. The remaining relations between £ c , £b and x^, s n _i are then given by energy and 
rapidity conservation requirements. The derived relations are explicitly listed in Appendix 

0- 

Combining all the derived rules for kinematic translation one finally obtains a trans- 



formation of Eq. 1.16: 



da^jAB -> X) 
dx a dxbdt n - \ds n -\ dip* 



(£c,6) daf^AB^X) 



J 



(s n -i,x b ) di a di b d^ 2 d^ c 



(2.3) 



where J 



is the Jacobian of the transformation derived in the Appendix [A|. 



(s„-i,x b ) 

An issue which deserves special consideration is the prescribed substitution t n —\ — 
m 2 = —/J- 2 , where the split particle virtuality fj? is the propagator virtuality shifted by 
the spectator mass m 2 . This prescription differs from the one of Collins ||, where the 
relation is directly t n ~\ = —fj, 2 and the spectator (propagator) mass shift is omitted. The 
reason for this modification is clear when one notes that the phase space limits [t~_ 1 ,t^_ 1 ] 
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of the t n -i parameter are functions of s n and invariant masses of the objects (particles) 



participating in a t-channel split of Equation 2.1 1£, |l^] and thus do not match the simple 



limits [m§, s n —i] of the virtuality fj? in the equation 1.16. Indeed, studies have shown that 
the presence of the cutoff 9(s n -i + i n -i) does not provide a sufficient solution since it only 
sets the upper integration limit to t^_ 1 — ► —s n -\ while the lower limit t~_ 1 can in certain 
instances be even smaller than m 2 . The reproduction of a logarithmic term of the collinear 
sine 



igularity In (^) (Eq. |TT|) with //| = /ig = s n _i: 



71 — 

/ 



g ( fl n-l+'n-l) rftH _ 1 = ln f Z£2zi ) ^ l n ( f^'l (2.4) 
*n-l V K-l I \ m c 



is thus not satisfied in when using t n -\ = —jj?. In order to resolve this issue one needs to 
return to the basics of the factorization procedure, where the actual collinear singularity 
(i.e. the corresponding logarithmic term) is isolated from the (integrated) hard process 
cross-section (for a nice example with massive particles see e.g. 0) and these logarithmic 
terms match with the required collinear logarithm ln (^i^j only in the high s n (hard 
center-of-mass) limit. The logarithmic collinear terms can subsequently be traced back to 
the propagator integral: 

dtn-l . ( 4-l~ m l \ , { S n -l\ 

^ = 7= ^2 sn^oo ln — T . (2.5) 



The expression of Eq. |2.5| is indeed found to match the logarithmic terms of |?]] exactly 2 . 
In order to reproduce the collinear cutoff one thus has to put — /j, 2 = t n -\ — m 2 which, 
combined with the factorization cutoff 0{pL 2 F — fi 2 ) = Q(s n -\ — fx 2 ), reproduces the required 
logarithm in the high s n limit: 



-(Ci-^D 



G ^ Sn K - - d/J. 2 = ln ( * n 1 — ] s n -» oo ln ( ) . (2.0) 

m§) 



2 Specifically the expressions of Eq. 17, page 14, whereby one has to write down the explicit expressions 
for the limits on t n -i and set the mass of the incoming particle corresponding to the incoming gluon to 
zero in order to reproduce the kinematical topology of the process studied in 
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2.2 Monte-Carlo Generation Steps 

A full Monte-Carlo event generation results 
using the following prescription: 

1. Generate the particle momenta and cor- 
responding phase space weight for the n- 
particle final state X and compute the full 
weight corresponding to the process of Eq. 



1.17 



2. Re-calculate the kinematic quantities of the 
n — > (n — 1) © 1 transformation as described 
above. 

3. Boost the whole system into the Collins- 
Soper frame g| of the n — 1 subsystem, 
i.e. the s n _i center of mass frame where 
the angle between the boosted hadron mo- 
menta P' A or —Pg and z'-axis now equals 
tan(±6>As) = l^xl/y^TI = \p£ \/y/s^Z{ (c.f. Fig 




Figure 3: The definition of the reference 
Collins-Soper frame as the rest frame of the 
subsystem s n _i with the transverse compo- 
nent k^_-y oriented in the x' — axis direc- 
tion. 



In other words, this transfor- 



2"AtS) — |'"n-ll/ V °n—l 

mation manifestly puts the transverse contribution due to the induced virtuality into 
the hadron momenta directed perpendicularly to the z'-axis. One can thus eliminate 
this virtuality by shifting the hadron momenta to the z'-axis while preserving the 
center-of mass energy s n -\. 



4. The Ti—l (hard) system corresponding to the Eq. |1.16| is then achieved by eliminating 
the particle c and boosting the remaining particles in the z-axis direction with the 
boost value of: 

fc(fc6* + ml) - tbiictbs + ml) 



13 idicibS + m 2 b ) + ibUdbS + m 2 c ) 
in order to restore the sub-system rapidity y. 



(2.7) 



5. Since boosts do not change the phase space weight the necessary modification consists 
of multiplying the phase space weight by: 



\{s n ,ml,ml)9(s n -i + t n -i)J- 



(2. 



(s n -i,x b ) 

and then obtaining the first subtraction weight by putting the reconstructed momenta 



into the Eq. 1.16 



6. An analogous procedure can be repeated to obtain the alternate kinematic config- 
uration (with the parton evolution of the other incoming particle) and the second 
subtraction weight. 

7. The final weight, after performing both subtractions is then passed to the event 
unweighting procedure. 
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As Chen, Collins et al. pointed out, H 10, 11, p^j , the described procedure of sub- 
traction is not equivalent to the standard subtraction schemes (e.g. MS) used to obtain 
the cross-sections for specific processes (see e.g. [23, and hence the data-fitted parton 



distribution functions with the corresponding evolution kernels (for example the widely 
used CTEQ5 PDF-s [^f]]). The relation between the procedure applied derived by Chen, 
Collins et al. for the massless case and (with inclusion of parton masses) applied above 
states that the correspondence between the MS scheme and the applied JCC scheme is 
given by relatively simple relations; an example for the expression of the quark i distribu- 
tion function involves a convolution of the gluon g distribution function and the g — > ii 
splitting kernel P g ^ii(z): 



zf i f I %,^) = zf^(z,^) (2.9) 

l 



a s (fl 2 ) f Z ms, 2 



„ ,Z. , / Z\ Z Z 



2, 



+ O (first-order quark terms) + 0(a 

These new distributions can in a reasonably straightforward manner be obtained by nu- 
merical integration. 

In order to complement the subtracted process one also has to generate a parton-shower 
evolved zero-th order process with (n — 1) particles participating in the hard process and 
an additional particle added by 'soft' evolution of the incoming particles. Since one is 
interested in the heavy initial state quarks this implies that one has to 'unresolve' one of 
the initial quarks back to a gluon, whereby an additional (anti) quark is added. The 
procedure to achieve this is straightforward |j [l(| (D], [0| an d complementary to the 
procedure described above, i.e. one has to perform the following steps in the Monte-Carlo 
algorithm: 

1. Generate the particles corresponding to the (n — 1) phase space topology, along with 
the momentum fractions £ c and of the incoming particles in the sense of light-cone 
components. Consequently, the invariant mass of the hard system is s n _i and the 
rapidity y is given by Equation 



Generate a virtuality /i 2 of the incoming heavy quark c, a longitudinal splitting 
fraction z for the branching of gluon a into the cc pair and an azimuthal angle (j) of 
the branching system. All the values are sampled from the Sudakov-type distribution: 

In case there are two quarks in the initial state that can evolve back to gluon and 
give the contribution of the same order (like e.g. bb — > Z° process) both virtualities 
are sampled and the quark with the higher one is chosen to evolve. 
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Subsequently, the four-momenta of the participating particles are reconstructed requiring 
that the subsystem invariant mass s n _i and the rapidity y are preserved; the construction 
is of course identical to the one used in the n — ► (n — 1) © 1 transformation given in Section 




2.1. A point to stress is that a kinematic limitation arises on the allowed x a and thus 
z = |^ values due to the requirement \p^\ 2 > 0. The latter condition gives the minimal 
value of the invariant mass of the n-particle system s n = x a XbS + m^ (taking into account 
that one of the incoming particles is a gluon, m a = 0) with: 

( W ) mm = { -^p^ ( (g - 1+ /i / f" + \i 1 : d + ^ i , (2.H) 

which combines with the rapidity y conservation requirement into: 

x 2 > (x a x b s) 2 min {(x a x b s) min - {m 2 c + 1?)) ^ 
a ~ s [(s n _i + fi 2 )e~ 2 y (x a x b s) min - m\{m\ + fi 2 )] 

In the massless limit the above expression translates back into the requirement x a > £ c 
or equivalently z < 1. In practice (i.e. Monte-Carlo generation) this thus means that a 
certain fraction of generated topologies have to be rejected and/or re-generated until the 
above conditions are met (or equivalently, that the z (or x a ) sampling limits have to be 
shifted). 
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3. Examples of the Algorithm Implementation 



Three examples of the procedure described in this paper have been developed: The asso- 
ciated Z°b production process and the 't-channel' and 'tW-channel' single top production 
processes, both expected to be observed at the LHC. The processes were implemented in 
the AcerMC Monte-Carlo generator jT?]]. Due to the subtraction terms a fraction of event 
candidates achieve negative sampling weights and unweighted events are produced with 
weight values of ±1 using the standard unweighing procedures. 

3.1 Associated Drell-Yan and b-quark Production 

The (Drell-Yan) lepton pair production associated with one or more heavy quarks repre- 
sents an important irreducible background component in the Higgs boson searches at the 
LHC. If the Higgs mass is around 130 GeV, then a promising decay channel is H — > ZZ* — > 
4£ where Z* represents an off shell Z. The production of a lepton pair with one or more 
heavy quarks is a background if the heavy quarks decay leptonically. Fig. ||| shows the 
relavent diagrams through order a s ; in each case there is a b and b in the event, at order 
a® both arise as fragements of the incoming beams. 




Figure 4: Representative Feynman diagrams for the Drell-Yan with associated b-quark production 
process for (from left to right): Order ai°\ 



order ai 1 ^ and order ai 1 ^ subtraction term. 



The cross-sections obtained for the leading order process bb — > Z — > ^ + fjT , next- 
to-leading order process gb — > Zb — > and the subtraction process (g — > bb)b — ► 
Zb — > in the LHC environment (proton-proton collisions at y/s = 14 TeV) are 
given in the Table [I], both for leading order PDF (CTEQ5L Q was used) and the PDF- 



s evolved according to the Collins prescription (c.f. Equation 2.1C , labeled JCC), along 
with the cross-section for the order as gg — ► Zbb — > fx + fi~bb process. Separate cross- 
section contributions for the next-to-leading order process and the subtraction term are 
given for convenience; in the Monte-Carlo procedure developed in this paper the events are 
generated according to the differential cross-section corresponding to the 'hard' order a!p 
process, i.e. the next-to-leading calculation with the subtraction terms subtracted on an 
event-by-event basis. 

The differential distributions of the virtuality fi and the transverse momentum {pt) 
distribution of the b-quark with the highest px of the two (produced either in the hard 
process or in the subsequent Sudakov showering) are shown in Figure ||, whereby also the 
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Mt or 1 £*<i£; 




cr CTEQ5L,/io=m z LP°J 


°JCC,^o=m z |P U J 


bb^Z^ 


g + g 


57.9 


39.9 


gb^Zb- 


-> g + g~b 


72.1 


60.0 


(g - bb)b 


— > Zb — > g + g~b 


73.3 


60.9 


£ 




56.7 


39.0 




— > g + g~bb 


22.8 


22.8 



Table 1: The process cross-sections for the leading order process bb — > Z — > g + g~ , next-to-leading 
order process <?6 — * Z6 — * g + g~b and the subtraction process (<? — * 66)6 — > Z6 — > g + g~b in the 
LHC environment (proton-proton collisions at = 14 TeV) are listed. The b-quark mass is set 
to mj = 4.8 GeV and the factorization and renormalization scales are set to the Z° invariant 
mass squared. In addition, the order ai 2 ^ gg — > Zbb — > g + g~bb process cross-section is shown for 



comparison. The cross-sections are given for the LO CTEQ5L 26 and the derived JCC PDFs. In 
the Monte-Carlo procedure the next-to-leading process weights are combined with the subtraction 
weights on the event-by-event basis as described in the text. 

/o\ _ _ 

Pt distribution of the b-quark with the highest px from the order Qs gg — > Zbb — > g + g,~bb 
process is shown. 




Figure 5: The differential distributions of the virtuality g and the transverse momentum (pt) 
distribution of the b-quark with the higher virtuality (produced either in the hard process or in the 
subsequent Sudakov showering) are shown for the calculations using the LO CTEQ5L parton density 
functions. The next-to-leading order process entries contain the subtraction terms as calculated on 
event-by event basis. 

In the results in Fig. |5| one can observe a smooth distribution of the virtuality g over 
the full kinematic range as the result of the implemented matching procedure. It is manifest 
that the cutoff on the b-quark virtuality g and the resulting subtraction contribution do 
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(2) 

not map to the pt distribution in a simple way. An interesting result is that the order 
gg — > Zbb — > fi + (i~bb process pt distribution of the b-quark seems to be quite close to the 
result of the merging procedure in the high kinematic range as one could indeed expect if 
the perturbative calculations are to be consistent in the perturbative regime. In the low 
Pt region the gg — > Zbb — ► fi + fi~bb process undershoots the expected distribution of the 
merged oi°^ © oi processes which is to be expected since in this case the non-perturbative 
contributions prevail. One has to keep in mind when comparing the two results results that 
in the derived calculation the other incoming b-quark is still effectively on-shell, i.e. its 
virtuality and branchings are obtained solely from the Sudakov showering, where as in the 

(2) 

as process both incoming b-quarks are treated as propagators in the full perturbative 
calculation. A further improvement would certainly be to repeat and extend the procedures 

(2) 

derived in this paper to include the full order a s calculation. 

From the results one can also see that the use of JCC evolved PDF-s significantly 
reduces the cross-section of the leading-order process with respect to the values obtained 
using the CTEQ5L ]2(| PDFs and to a lesser extent the cross-sections of the next-to-leading 
and the subtraction contribution, since the latter two include only one b-quark and one 
gluon in the initial state and are thus less affected by the change in the b-quark PDF 
evolution. 



3.2 The 't-channel' Single Top Production Process 

The 't-channel' single top production mechanism is of importance at the LHC since it 
provides a clean signal for top quark and W boson polarization studies. The final state 
consist of a t, W, and b as illustrated in Fig. |6| 





Figure 6: Representative Feynman diagrams for the 't-channel' single top production process for 
(from left to right): Order ai°\ order Qs 1 ' and order a^ 1 ' subtraction term. 



The cross-sections obtained for the leading order process bq — ► tq' , next-to-leading or- 
der process gq — > tq'b and the subtraction process (g — > bb)q — > tq'b (and charge conjugates) 
in the LHC environment (proton-proton collisions at y/s = 14 TeV) are given in the Table 
|2| both for leading order PDF (CTEQ5L |2(| was used) and the PDF-s evolved according 
to the Collins prescription (c.f. Equation |2.10 , labeled JCC) for the scale choices = m t 
and fiQ = 60 GeV. The differential distributions of the virtuality \i and the transverse 
momentum (pr) distribution of the b-quark (produced either in the hard process or in the 
subsequent Sudakov showering) are shown in Figures |7| and H Separate cross-section con- 
tributions are given for convenience; in the Monte-Carlo procedure developed in this paper 
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the gq — > tq'b events are generated according to the differential cross-section corresponding 
to the 'hard' order process, i.e. the next-to-leading calculation with the subtraction 
terms subtracted on an event-by-event basis. 



Process 


a CTEQ5L, [pb] 
M0= m t 


<t jcc, [pb] 

M0= m t 


a CTEQ5L, [pb] 
/i =60 GeV 


a jcc, [pb] 

MO =60 GcV 


bq -> tq' 


222.2 


187.8 


178.1 


138.7 


gq -> tq'b 


156.2 


154.2 


188.2 


184.4 


(g - bb)q - 


-» i<?'6 140.1 


138.2 


102.8 


100.5 


£ 


238.3 


203.8 


263.5 


222.6 



Table 2: The process cross-sections for the leading order process bq — > £(?', next-to-leading order 
process — ► and the subtraction process (g — > 66)g — > (including the charge conjugates) 
in the LHC environment (proton-proton collisions at i/s = 14 TeV) are listed. These inclusive 
cross-sections include all top (and W^) decay channels, the b-quark mass is set to m& = 4.8 GeV 
and top-quark mass to m t = 175 GeV with the factorization and renormalization scales set to the 
top mass values no = m t and no = 60 GeV. The cross-sections are given for the LO CTEQ5L and 
JCC PDFs. In the Monte-Carlo procedure the next-to-leading process weights are combined with 
the subtraction weights on the event-by-event basis as described in the text. 

From the results in Fig. |j] and Fig. || one can observe that the applied procedure 
produces a very good match of the processes in the combined distribution of the b-quark 
virtuality fx resulting in a smooth (almost seamless) transition in the vicinity of the cutoff. 
As one can also observe the cutoff on the b-quark virtuality /i and the resulting subtraction 
contribution do not map to the pt distribution in a trivial manner; hence one can surmise 
that the simple methods involving adding of the processes based on pp distribution cuts 
probably give erroneous predictions. 

This procedure can in unmodified form be applied to the full 2^4 and 2^4 matrix 
elements bq — ► tq' — ► Wbq' — > ff'bq and gq — > tq'b — > Wbq'b — > ff'bqb including top quark 
decays and has as such also been implemented in the AcerMC Monte-Carlo generator. 
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Figure 7: The differential distributions of the virtuality [i and the transverse momentum (pt) 
distribution of the b-quark (produced either in the hard process or in the subsequent Sudakov 
showering) are shown for the calculations using the LO CTEQ5L PDFs and the showering (factor- 
ization) scale set to = 60 GeV. The next-to-leading order process entries contain the subtraction 
terms as calculated on event-by event basis. 




Figure 8: The differential distributions of the virtuality /j, and the transverse momentum (pr) 
distribution of the b-quark (produced either in the hard process or in the subsequent Sudakov 
showering) are shown for the calculations using the JCC PDFs and the showering (factorization) 
scale set to ^ = m t- The next-to-lcading order process entries contain the subtraction terms as 
calculated on event-by event basis. 
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One of the interesting results is that the use of JCC evolved PDF-s significantly reduces 
the cross-section of the bq — ► tq' process with respect to the values obtained using the 
CTEQ5L PDFs and to a lesser extent the cross-sections of the gq — * tq'b and the subtraction 
contribution, since the latter two contain only the light quarks and gluons which are less 
(or in case of gluons not at all) affected by the change in PDF evolution compared to the 
b-quark PDF. 

A somewhat cruder method of merging different order processes has for the 't-channeF 



single top production already been implemented a while ago in the program ONETOP [27]; 
the results from the two methods are compatible within the differences of the methods used 
in both implementations. It is instructive to understand where the differences between the 
two procedures originate. Specifically, in ONETOP the subtraction term is higher than 
the cross-section of the (2 — > 3) gq - 
cross-section of the process (g — > bb) 



tq'b process; from the Fig. 3.6 in 27] the subtraction 
— > tq'b is of the order of about 190 pb compared to the 
140 pb one gets using the procedures described in this paper. The difference originates in 
part in the massless approximation of the participating particles implemented in ONETOP 
and can be traced back to the fact that the subtraction term in ONETOP (as given in 
Appendix D in p7| ) is calculated from the integrated parton density function correction 
(the first-order term in Eq. 1.12 of this paper) coupled to the zero-th order bq — > tq' cross- 



section, since the virtuality is already integrated over into the log-^j. In contrast, in the 
present work the procedure is more complex and requires identifying the virtuality of the 
2^3 process in the massive calculation. In addition, in ONETOP the spectator energy 
fraction is kept unchanged whereas in the new procedure only the rapidity constraint is 
used instead. The ONETOP calculation using the integrated PDF correction has been 
repeated as a check and gives a subtraction term with the value of about 185 pb which is 
consistent with the ONETOP results. 

3.3 The 'tW channel' Single Top Production Process 

The 'tW-channel' single top production mechanism is also of importance at the LHC since 
it provides a clean signal for top quark and W boson polarization studies. The process is 
illustrated in Fig. 0. 






Figure 9: Representative Feynman diagrams for the 'tW-channel' single top production process 
for (from left to right): Order a$ , order ai 1 ^ and order c4 subtraction term. 

The cross-sections obtained for the leading order process gb — > t(W —►)//', next-to- 
leading order process gg — ► t(W —*)ff'b and the subtraction process (g — > bb)g — ► t(W — ► 
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)ff'b (and charge conjugates) at the LHC. are given in the Table ||, both for leading order 
PDF (CTEQ5L (2f| was used) and the PDF-s evolved according to the Collins prescription 
(c.f. Equation p. 10 , labeled JCC) for the scale choice fio = 60 GeV. The differential 
distributions of the virtuality fi and the transverse momentum {pt) distribution of the 
b-quark (produced either in the hard process or in the subsequent Sudakov showering) are 
shown in the Figure . Separate cross-section contributions are given for convenience; 
in the Monte-Carlo procedure developed in this paper the gg — > t(W —►)//'& events are 
generated according to the differential cross-section corresponding to the 'hard' order ctg 
process, i.e. the next-to-leading calculation with the subtraction terms subtracted on an 
event-by-event basis. Note also that in this case there are indeed two subtraction terms, 
one for each incoming gluon. 

From the results in Fig. [R^ one can again observe that the applied procedure produces 
a very good match of the processes in the combined distribution of the b-quark virtuality 
H resulting in a smooth (almost seamless) transition in the vicinity of the cutoff. As 
one can also observe the cutoff on the b-quark virtuality [i and the resulting subtraction 
contribution do again not map to the pt distribution in a trivial manner; again therefore 
one expects that the simple gluing methods of the processes based on px distribution cuts 
most probably give erroneous predictions. 



Process 




O"CTEQ5L,/i =60 GeV [pb] 


CJCC>o=60 GeV [pb] 


gb -► t(W - 


-OA"*!* 


5.9 


4.7 


gg - t(W - 


->)fjbv tl b 


6.0 


6.0 


(g - bb)g - 


-> t(W ->)/iP M & 


3.1 


3.1 


E 




8.8 


7.6 



Table 3: The process cross-sections for the leading order process gb — ► t(W —►)//' , next-to-leading 
order process gg — > t(W —>)//'& and the subtraction process (g — ► bb)g — ► t(W —>)ff'b (including 
the charge conjugates) in the LHC environment (proton-proton collisions at y/s — 14 TeV) are 
listed. These inclusive cross-sections include all top decay channels, the associated 1T ± decays into 
muon and neutrino and the b-quark mass is set to rat = 4.8 GeV and top-quark mass to m t = 175 
GeV with the factorization and renormalization scales set to the /xq = 60 GeV. The cross-sections 
are given for the LO CTEQ5L and JCC PDFs. In the Monte-Carlo procedure the next-to-leading 
process weights are combined with the subtraction weights on the event-by-event basis as described 
in the text. 

The diagrams of the process gg — * tWb — > WWbb — > ffffbb are in fact just a subset 
of 31 Feynman diagrams which have a WWbb intermediate state (and which also include 
the tt production). Accordingly, the derived subtraction procedure has in, AcerMC, been 
applied to the processes having the full set of Feynman diagrams and the above plots and 
values should be considered only as the validation of the procedure in case of the 't W 
channel' single top production. 
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Figure 10: The differential distributions of the virtuality /i and the transverse momentum {pr) 
distribution of the b-quark (produced either in the hard process or in the subsequent Sudakov show- 
ering) are shown for the calculations using the LO CTEQ5L PDFs and the showering (factorization) 
scale set to /io — 60 GeV. The next-to-leading order process entries contain the subtraction terms 
as calculated on event-by event basis. 



4. Conclusion 



It has been demonstrated explicitly how to deal with the case where a particle of interest 
can be produced from a partonic hard scattering process or as a remnant of an incoming 
hadron beam. Examples of the proceedure have been provided and contrasted with the 
more ad-hoc proceedures used previously to prevent double counting. 
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A. Kinematic Relations 



The requirement (2) in the list of Subsection 2T gives the equivalence £ a = x a- The 
remaining relations between £ c ,£b and Xb,s n -i are then given by energy and rapidity 
conservation requirements and are thus given by the conditions: 

2 2 

Sn-i = (p c + Pb? = m 2 c + m 2 b + 2(p+p- + p-p+) =m 2 c +m 2 b + £ c £ b s + (A.l) 
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where p c 
and: 



with: 



and: 



(p+ ffTp-) = ^cPt,^,0r) and p b = (p+ T ,p~) 

1 - / 6c6f>S + 



In 



"n-l 



v n-l 



111 



66 



"I - ~ In i 

2 V £ c &>s + m% 



v n-l 



1 + /9 

1-/3 



\l-l 



v n-l 



1 + /A _ x a (x a XbS + ml) 



1-/3/ rE fc (rE a x 6 s + m2) 



(A.2) 



(A.3) 



v n-l 



1 

71 



g n + s n -i - mj 



t n _i - mg - s n _i 



+ 2 



(A.4) 

which can be inverted to give the expressions for £ c and as functions of s n , x^, ■ ■ ■ A 
further simplification in derivation can be achieved by introducing another set of variables 
f = £ c • £{, and y = l/2m(£ c /£f,) with the Jacobian of the transformation J ^§4^ = 1 and 



subsequently: 



Y s { - ( m c + m b)) + y/K s n-i, 



1 in ^ 

2 Un-1 



rs + mi 



(A.5) 
(A.6) 



As one can observe the r is only a function of s n _i so the only remaining term to compute 



in the Jacobian of the transformation 



J xi^k) is d y/ dXb: 



j 



(i~,y) 



(s n -i,x b ) 



df 



ds 



n-l 



dy 



dxh 



J~(s n , s n —\ , in— i j 3Ca 1 m a , mb, m c ) 



(A.7) 



with T being a lengthy function of the listed parameters and therefore omitted. In the 
massless approximation the above expression reduces to: 



J 



(6c, 6ft 



(s n -lj x b) m^O 



2(x a XbS + t 



n—lj 



2XbS 



(Ai 



which is in agreement with the expressions derived by Chen, Collins et al. |)[ [T(], 11, 12] 
Analogous expressions can trivially be obtained also for the split of the other parton. 
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